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ABSTRACT 

The kinetics of irreversible annihilation of charged particles perform- 
ing overdamped motion induced by long-range interaction force, F(r) ~ 
r _A , is investigated. The system exhibits rich kinetic behaviors depend- 
ing on the force exponent A. In one dimension we find that the densities 
decay as t~ l ^ 2+x) and t~ l ^ 1+2X) when A > 1 and 1/2 < A < 1, respec- 
tively, with logarithmic correction at A = 1. For A < 1/2, the asymptotic 
behavior is shown to be dependent on system size. 
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1. INTRODUCTION 



The kinetics of two-species diffusion-controlled annihilation reaction, A + B — » 0, be- 
tween uncharged particles has been a subject of extensive research for almost 20 years [1] . 
For sufficiently low spatial dimension, d < 4, even under homogeneous initial conditions, 
large-scale heterogeneities arise that invalidate classical kinetic laws. Much less is known 
about annihilation reaction between charged particles with long-range power-law interac- 
tion, F(r) ~ r~ x . An important case of Coulomb interaction (A = d — 1 in d dimensions) 
has been treated in a few studies [2-5] for d = 2 and 3. However, some of these works were 
based on unjustified approximations while others were based solely on numerical simula- 
tions, so their results are also uncertain. (Note that even for annihilation of uncharged 
particles in three dimensions the asymptotic regime is hardly reached on modern comput- 
ers). We, therefore, see that the Coulomb case still deserves further investigation. Other 
values of the interaction exponent A also naturally appear in applications with particles 
being dipoles, defects, vortices, monopoles, disclinations, etc. One important example is 
the quench of a one-dimensional Ising system from a disordered state to an ordered state. 
If spins interact via long-range potential [6], the Hamiltonian may be expressed in terms 
of interacting domain walls [7]. There are two types of domain walls in the system, the 
domain walls with "up" spins to the right and "down" to the left (A- walls), and vice 
versa (S-walls). Thus, an alternating domain wall sequence ...ABABAB... is formed. 
Domain walls annihilate upon colliding, A + B — > 0, but since the alternating structure 
persists in time, the reaction process is, in fact, equivalent to the single-species annihila- 
tion, D + D — > 0. This system has been recently investigated [7-9], and it was shown that 
particle concentration decays as t _1 /( 1+A ). 

In this paper we consider a truly two-species annihilation model where the initial 
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distribution of interacting particles ("charges") is random (Poissonian). The forces between 
charges are assumed to be proportional to r _A , with similar charges repelling each other 
and dissimilar attracting each other. Compared to the single-species case, the two-species 
annihilation exhibits more rich kinetic behavior including the dependence on the system 
size. In this study we focus on one-dimensional systems which allow us to find rather 
convincing numerical support for our scaling predictions. 

The paper has the following structure: In the next Section we formally introduce the 
model, an ensemble of interacting particles in one dimension with overdamped dynamics. 
Then, relying on heuristic arguments, we obtain density decay exponents for different values 
of A. In Section III, we present the results of numerical simulations. Finally, in Section 
IV we discuss possible generalizations of the model including higher dimensionality and 
ballistic motion, and make general conclusions. 



We consider two-species systems containing A- and B-type particles with charges +1 
and —1 for "particles" and " antiparticles" , respectively. Particles of both species move 
continuously in one dimension and interact via long-range force, F = qq' jr x (for charges 
q and q' separated by the distance r). Initially, A- and S-type particles are randomly 
distributed with equal concentrations, for simplicity we put them equal to 1. 

Total force acting on the i th particle is equal to a sum of pairwise forces: 



where qu = ±1 and Xk are charge and coordinate of the k particle. We will ignore 
particle inertia; i. e., motion of particles is assumed to be overdamped. Therefore, the 



2. THE MODEL AND SCALING ARGUMENTS 
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velocity of each particle is proportional to the total force Fj acting on it, Vi = fxFi. (In 
the following, we set the mobility [i equal to 1). We will also ignore particle diffusion; that 
is, we will assume that the drift dominates the random walk effects. When two dissimilar 
particles collide, both of them irreversibly disappear; collisions between particles of the 
same species are impossible because of repulsion. To summarize, we consider two-species 
annihilation of particles undergoing overdamped noiseless motion. We will see that for a 
sufficiently large force exponent (A > 4 in one dimension), the noise actually dominates 
the drift, and thus well-known diffusion-controlled kinetic behavior [1] emerges. However, 
for small A, we expect that the long-time behavior is correctly described by our noiseless 
model. We will also briefly discuss a model where the motion is ballistic (Sec.V). 

Let us now consider time evolution of the system. We cannot a priori expect that a 
mean-field description holds in low dimensions, especially in one dimension. Remember 
that the breakdown of the mean-field behavior in reaction-diffusion models is generally 
attributed to the formation of single-species domains [1] . We assume that the same takes 
place in our model (at least in one dimension, the formation of domains is inevitable). 

Suppose that at time t the length of a typical single-species domain is L(t). It means 
that an average number of particles in such a domain is equal to initial imbalance of 
majority and minority species on the length L(t), which for Poissonian initial distribution 
with the density one is of the order of y/L(t). Therefore, the concentration n(t) in a typical 
domain behaves as 

n(t) ~ l/v^L^t). (2) 

To get an insight on how a typical domain length changes in time, we consider motion 
of a single particle A on a domain edge (Fig. 1). Here we assume that the system is entirely 
formed of well-defined domains of typical length L. The total force acting on A from the 
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particles on its left may be evaluated as: 



M ^ M 
< ™A 



M 



(3) 



Here M is a typical number of particles in a domain, M ~ \/L. Each sum in the left-hand 
side of Eq. (3) expresses a contribution to the net force from a particular domain; to get 
the total force these contributions have been added. It is clear that the force exerted on A 
from all the particles to the right may be calculated in exactly the same way. To simplify 
the matter even further we assume that Xj ~ R x j, where R is an average interparticle 
separation; R = L/M. Rewriting Eq. (3) through R and M gives 



F ~ — 
R x 



M M M 



i 



(4) 



Depending on the value of the force exponent A, different situations appear. For A > 1, 
the first sum converges to a finite value for M — > oo while the other sums approach to 
zero as M~^~ x ' which means that only charges from the left and right nearest neighbor 
domains essentially contribute to the total force. For < A < 1, all sums diverge as 
M 1 ~ x as M — > oo. The total force, being the sum of monotonically decreasing sign- 
alternating terms, is of order of the contribution of the first domain. In the borderline case 
A = 1, the first sum diverges as logM while the others terms are monotonically decreasing, 
sign-alternating, and finite. The dominate contribution is again provided by the nearest 
domains. To sum up, 



F ~ — r- X < 

R x 



( M x ~ x if 0<A<1, 
logM if A = l, (5) 
1 if A > 1. 

On the other hand, a typical rate of change of the domain length is of order of the 
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velocity of any of its edges. Recalling that R ~ M ~ V% we obtain 

' L i/2-\ if <A<1, 

dL/dt~F~ \ L-V^ogL if A = l, (6) 
> L" A / 2 if A > 1. 

Solving (6) for L{t) and using relation (2), we finally write for the density decay 
asymptotics: 

I.* if 0<A<1, 
(Hog*) - * if A = l, (7) 
if a>1. 

These results are, in fact, correct only for 1/2 < A < 2. The upper bound follows from 
comparison of the random walk length, Lrw ~ t 1 / 2 , with the drift length, L ~ t 2 /( 2 + A ) 
when A > 1. For A > 2, Lrw 3> L, so a pair of charges can escape annihilation through 
a random walk, and therefore the diffusion controls the dynamics. Thus, for A > 2, the 
diffusion-controlled asymptotic behavior, n ~ t -1 / 4 , is expected. The lower bound stems 
from the fact that an average force acting on any particle in the infinite-particles ID 
system becomes infinite for A < 1/2. It can be shown rigorously by deriving a Holtsmark- 
like [10] force distribution (we leave it for Appendix A). Here, we provide more qualitative 
arguments which take into account the finiteness of the system. First, we note that in 
calculation of the total force in (3), we implicitly assumed that the system is perfectly 
ordered — it consists of similar domains of A and B particles; i. e., the total charge of 
the first domain is equal, up to the sign, to the total charge of the second domain, etc. In 
particular, it means that for a system depicted in Fig. 1, the overall charge to the left of 
the test particle A is —1, and the overall charge to the right is zero. However, this picture 
is a "mean-field" in spirit; hence, it can lead to erroneous results for truly random systems. 
Fluctuations in initial charge distribution in the system with iV particles produce the net 
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charge of the order to the left and to the right of the test particle. Qualitatively, we 
can estimate the effect of charge imbalance by putting v^V equidistant charges of one sign 
to the right and the same amount of the opposite charges to the left. (When we consider 
finite systems of charges we implicitly assume that they satisfy the neutrality condition; 
generally, the total net charge determines the long-time behavior). The initial size of the 
system is N, so the distance between nearest charges left in the system is \/~N ; thus, the 
force F N due to the charge imbalance is 

i=i J 

While Fn — > as N — > oo, this force does not affect the dynamics of the model. Thus, 
for A > 1/2 our previous estimate, F ~ L 1//2_A , gives the dominate contribution, and 
the size-independent dynamics (7) emerges. However, for A < 1/2 the total force acting 
on a particle grows with system size even for overall neutral systems; therefore, it should 
dominate over the "regular" force (3) and control the dynamics of the system. It seems 
reasonable to assume that at the early stages of time evolution, when the distribution of 
the particles is still almost random, the motion of domain edges is controlled by the force 
(8). Repeating the steps used in deriving Eq. (7), with F ~ jV 1 / 2_A instead of Eq. (6), we 
obtain for the density decay 

n(t) ~ N^t~^. (9) 

However, this estimate may become inapplicable on the later stages of evolution. Indeed, 
the dynamics described by Eq. (9) is extremely fast since it is size-dependent, so the charge 
distribution that emerges can be significantly different from the Poissonian; as a result, 
our assumption about the type of randomness of the particle distribution could become 
less and less appropriate. 
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3. SIMULATION RESULTS 



To check our heuristic predictions, we have performed numerical simulations for A = 
1, 0.75, 0.5, 0.25, and (the Coulomb case A = turns out to be special; it will be discussed 
separately in Appendix A). Our system initially consisted of 10000 particles of each species 
randomly distributed with concentration 1. First, the net force (1) is calculated for each 
particle. We compute all the forces directly without applying any multipole-like expansion 
that could be useful in many dimensions [12]. With the particle velocity equal to the total 
force, we employ a simple Euler update procedure for each time step: Axi = FiAt. The 
selection of time interval At was merely experimental. Since on the last stages of evolution 
simulations run very fast (few particles are left), we, unlike [7-8], keep At constant during 
a run. Finally the results are averaged over 10 runs. The selection of boundary conditions 
does not seem to affect the results of simulations of the two-component system with overall 
neutrality, except for, maybe, its latest stages. We ran the simulations for A = 0.25 with 
both periodic and open boundary conditions; the results for concentration coincided within 
a statistical error. Hence, we use open boundary conditions for all further simulations. 

In Fig. 2 we plot, for various A, the concentration of either species vs. time. As 
we anticipated, for 1/2 < A our predictions for the decay of concentration are in good 
agreement with the results of simulation. Indeed, the average slope of log(n(£)) vs. logt 
plot for A = 1 is —0.322 [the heuristic argument gives (tlogt) 1 ], and —0.401 (compared 
to —2/5) for A = 3/4. Performing numerical simulations with a twice smaller system (5000 
pairs of particles) , we did not find any significant system size dependence for these values 
of A. As for A < 1/2, no power-law behavior can be observed for density decay. However, 
at the early stages of evolution, local exponents are close to —1/2 as it follows from Eq. (9); 
as time goes on, the density decay rate increases. 
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4. DISCUSSION 



Considering the dynamics of domain interfaces and using simple heuristic arguments, 
we have predicted the asymptotic density decay in a two-species annihilation reaction 
system with long-range interaction and overdamped motion. This approach can be gen- 
eralized on similar systems in an arbitrary spatial dimension d. Following the line of 
reasoning described in Sec. II, one can obtain: 

' t -2=j+2x jf d/2 < A < d, 

(tlogt)"^ if X = d, ( 10 ) 

W~2+X jf \> d. 

For A < d/2, we expect size-dependent kinetics. For A > 2, as in the Id case, let us 
compare diffusion and drift length scales, L RW ~ t 1 / 2 and L ~ £ 2 /( 2 + A ) ? and conclude that 
the random walk dominates the drift; hence, the diffusion-controlled behavior, n ~ t _d / 4 , 
is expected. Thus, the regime described by the lower line of Eq. (10) does not appear for 
d > 2. Note that for truly Coulomb systems, A = d — 1, the scaling prediction of Eq. (10) is 
n ~ t -1 ; i. e., the classical kinetic law. This is expected to arise when d/2 < A = d — 1 < 2, 
i. e., 2 < d < 3. For d < 2, size-dependent kinetics is anticipated; while for d > 3, the 
Coulomb interaction becomes irrelevant and the diffusion-controlled behavior emerges. 
Another interesting example is the system of D-dimensional Coulomb charges confined 
to the hypersurface; i. e., A = d = D — 1. In this case, Eq. (10) predicts logarithmic 
corrections to the power law behavior. 

However, it is not clear whether the concept of unpenetrable (untransparent) domains 
with continuous boundaries is still applicable for d > 1. Competition between screening, 
which makes long-range interaction effectively short-range, and annihilation may also sig- 
nificantly affect the behavior of the system. We attempted to simulate a 2d system on a 



lattice. Either because of the large effective diffusion, which is inevitably introduced by the 
discrete nature of the lattice model, or for some deeper reason, we were unable to observe 
any scaling-like behavior. Since many-dimensional continuous many-body simulations are 
still computationally challenging, we leave this problem for the future. 

More can probably be done with one-dimensional systems as well. One can try to find 
exact results for some specific values of the force exponent A. For the Coulomb system, 
A = 0, we indeed succeeded in finding some properties analytically (see Appendix A and 
Ref. [11]), but we still could not find a complete solution. Another extreme case, A — > oo, 
is also theoretically challenging. (In fact, the diffusion determines the dynamics for A > 4 
so one should consider a strictly noiseless system). The dynamics in the A — > oo limit 
is extremal: One picks the pair of nearest neighbors that are closest to each other and 
removes it if the charges are dissimilar, or recedes them if the charges are similar. The 
receding is stopped when the distance reaches the second minimal intercharge distance. 
Then, if this second closest pair contains dissimilar charges, it is removed while the first 
pair continues receding; if the second pair is also the same-species, both pairs recede. If the 
initial sequence is alternating as it takes place with domain walls in the quench process, 
one always removes the closest pairs, and the domain-size distribution function approaches 
the scaling form. This model turns out to be completely solvable [13-15,8]. It would be 
very interesting to study a more complex version of the extremal dynamics that arises from 
the present two-species annihilation process. 

Finally, we discuss a model where motion of the particles is ballistic; i. e., it is de- 
scribed by Newton's laws. It is clear enough that the only change one needs to make in the 
above approach is to put d 2 L/dt 2 instead of dL/dt into the left hand side of (6). Assuming 
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that initial velocities are irrelevant, we obtain for the concentration 




1 



if 



A < 1/2, 



2 



1/2 < A < 1, 



t 1 + 2X 



if 



n(t) ~ < 



(tlogt) 



2 
3 



if 



A = l, 



(11) 



2 




t 2 + A 



if 
if 



1< A < 2, 
A > 2. 



When A > 2, the inertia dominates the drift and, therefore, the ballistic-controlled asymp- 
totic behavior [16], given by the last line of Eq. (11), follows. 

Numerical simulations performed for A = 0.75 showed that the concentration decays 
as t -0 - 79 compared to t~ 4 / 5 as it follows from Eq. (11). 
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APPENDIX A: A = 



The A = case, corresponding to the "real" Id Coulomb force has several peculiar 
features. Since the forces between particles do not depend on the distance, and particles 
disappear in pairs, the net force acting on a particle is constant throughout its life. It 
means, that velocity of any particle is constant and equal to the difference between total 
charge to the left and total charge to the right of it multiplied by a charge of the parti- 
cle. Taking into account that initial distribution of particles is Poissonian, one can easily 
describe the behavior of the system if the probability distribution for " charge imbalance" 
were known. To get an insight about this charge imbalance distribution, we look at our 
configuration of charges as on a Id random walk (RW) (Fig. 3). Step up corresponds to 
a positive charge, step down - to the negative; since the system is overall neutral, the 
RW returns to the origin after 2N (the size of the system) steps. The net force acting on 
a particle is equal to the "height" h, positive or negative, of the corresponding point in 
the RW picture. The joint distribution function W2N(h,2L) tells us how many segments 
(loops in the RW terminology) of the length 2L, starting and ending at the "height" h 
from the origin, exist in a RW coming to the origin (not necessarily for the first time) after 
2N steps. Knowing this function, one can readily calculate a life expectation time for each 
particle, and therefore the concentration decay rate: 



Here L is initial system length and Pj(x) = x^e~ x / j\ the Poisson distribution function. So 
far we have been able to determine another function, W / 2jv(2£), which gives the probability 
that the maximum length of "zero-height" segment in a system with 2N charges is 2L 
[11]. This maximum length segment determines the lifetime of the whole system, which is 




N N-L 
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shown to be proportional to its size 2N. Moreover, the life time distribution function has 
a remarkably rich structure (an infinite set of singularities, etc.] see [11]). 

Although we were not able to find the exact expression for the decay of concentration 
in this case, numerical simulations of this problem prove to be very simple. Instead of 
running a molecular dynamics algorithm, it is sufficient to calculate all the net forces once 
and find an annihilation partner for each particle; after that, we know the lifetimes for all 
the particles in the system. It enabled us to check an accuracy of our molecular dynamics 
simulation for our standard (10 4 pairs) system size and also to study the systems with up 
to 10 5 particles of each species. Even for these relatively big systems, we were unable to 
find any power-law behavior; the function that fits best our simulation results looks like 
n(t) ~exp(-0.051n 4 t). 



APPENDIX B: 

TYPICAL FORCE IN A SYSTEM OF CHARGED PARTICLES 

A problem of distribution of a typical force in a system of particles with Coulomb 
interaction was first studied by Holtsmark [10]. Using the same approach, we will show, 
that for arbitrary dimensionality and r~ x interaction, the typical (or mean) force is finite in 
an infinite system only for some range of the force constants A, specifically for A > d/2. We 
start from an expression for force distribution function W(F), which gives the probability 
that the force, acting on a "test particle" which we will put at the origin, is equal to F: 

N 

W(F) = (5(F-J2m)))av, (Bl) 

where fifj) = rj/r x+1 is the force exerted by a j th particle on the test one. Since we 
assume that spatial distribution of particles is random and independent, it is sufficient to 
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consider one-component system. Replacing 5-function by auxiliary integration over dk we 
obtain: 



W(F) = 



1 



(2tt)< 



exp (ik-F)S{k)dk 



(B2) 



with 



S(k) 



N 



N 



exp -^ik-fft) 11 



V \V 



1 



- ik ' f df 



N 



(S3) 



j = l / 3=1 

Here j = f/r x+1 , N is the number of particles in the system, and V is the volume. 
Rewriting (B3) in the form 



S(k) 



1 



(S4) 



and taking the thermodynamic limit, iV — > oo and F — > oo where A^/F = n is kept fixed, 
yields 

S(jfe) = exp ^-n J\l- e -^-7- A+1 j ^ . ( B5 ) 

After expanding the exponent in the integrand for large r and performing angular 
integration (which eliminates all odd-order terms) one finds that for system size R — > 
oo, the integral in the right-hand side of Eq. (B5) converges only if A > d/2 while for 
A < d/2 the integral diverges with the system size. More precisely, a considerable but 
straightforward computation yields 



S(k) = S(k) = exp -ntt d A(d, X)k 



(B6) 



for A > d/2. In Eq. (B6), fi d = 2ir d / 2 /T(d/2) denotes the surface area of the unit sphere 
in d dimensions and A(d, A) is the shorthand notation for the integral 



A(d, A) 



1 dz 
o Xzi +1 



i -r 



J d 

2 
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with T(z) being the Euler gamma function, and J u (z) the Bessel function. For A < d/2 
the integral in the left-hand side of Eq. (B5) grows as g^gzf^ R d ~ 2X k 2 . Returning back to 
Eq. (B2) we see that in the case of (A < d/2) the net force is given by 

w & - Tk*l e ^ d"- p - w^f 2 ^) dt <B7) 

In the limit R — > oo, we compute the integral in Eq. (B7) asymptotically to find 

W(F) ~ [(d - 2\)R 2X ~ d F] d ~ X exp [-const (d - 2\)R 2X ~ d F 2 ] . (B8) 

Eqs. (B7) and (B8) are valid for A < d/2. For A = d/2, R d ~ 2X /(d - 2A) should be 
replaced by logR. For d = 1 and R oc N, the net force has the Gaussian distribution, 
W(F) oc (l-2A) 1/2 A^ A " 1 /2 exp _ 2A)A^ 2A_1 F 2 ] , and therefore the typical force grows 
with size as N l / 2 ~ x in agreement with qualitative results of Sec. III. Note also that for 
A = 1/2 the typical force still grow with system size, F ~ \/\og N, and hence the density 
decay of the form n(t) ~ (log iV) -1 / 4 ^ -1 / 2 is expected. 
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Figure Captions 

Figure 1. Schematic illustration of domain structure. Typical domain has length L and consists 
of M = 6 particles, average distance between particles is R = L/N. 

Figure 2. Plot of concentration n(t) vs. time on a double logarithmic scale for A = l(o), A = 
3/4(A), A = 1/2( V )A = 1/4(0), and A = 0(*). 

Figure 3. RW representation of a two species neutral system. Number of particles of each 
species N = 17. A neutral segment of the length 2L = 16 is shown having h = 2 
uncompensated charges to the left and to the right of it. 
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